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RECENT DEVELOPMENTS IN FINITE ELEMENT ANALYSIS FOR TRANSONIC AIRFOILS* 

N. M. Hafez and E. M. Muraan 
Flow Research Co., Kent, Washington 

INTRODUCTION 

The prediction of aerodynamic forces in the transonic regiine generally 
requires a flow field calculation to solve the governing non-llnt^ar mixed 
elliptic-hyperbolic partial differential equations. Finite difference tech- 
niques have been developed to the point that design and analysis application 
*^®**t*ne, anc. continual Improvements are being made by various research 
groups. The principal limitation in extending finite difference methods to 
complex three-dimensional geometries is the construction of a suitable mesh 
system. Finite element techniques are attractive since their application to 
other problems have permitted irregular mesh elements to be employed. The 
purpose of this paper is to review the recent developments in the application 
of finite element methods to transonic flow problems and to report some recent 

results of our own study. In most cases, the reader is referred to the original 
paper for the details. ui.isxii«x 


Finite element methods have been quite sucessful when applied to elliptic 
problems, particularly in elasticity and structures. A straightforward appli- 
cation of these techniques to the transonic problem involving mixed elllptlc- 
hyperbollc equations with embedded shock waves has not proven successful. In 
general, either the finite element method must be modified to treat the tran- 
sonic flow equations or the transonic flow equations must be modified to fit 
the finite element method. For the latter approach, additional terms (with 
hopefully small coefficients) are added to the equations representing arti- 
ficial viscosity, artificial density or artificial time. The general mesh 
shapes used for elliptic problems result in matrices which are not well 
ordered. Thus either a direct matrix Inversion or an explicit Iterative method 
is required. Both of these approaches are computationally inefficient for 
realistic problems compared to the implicit methods which have been developed. 
Some work reported later in this paper linay relax this conclusion. Finally, the 
higher order shape functions used in elliptic problems do not apperr as attrac- 
tiv 3 for mixed problems with shocks. In this paper, we review the subject by 
generally grouping together methods follovd.ng the same general formulation. 

SYMBOLS 


.All symbols are dlmenslonals. 

A matrix operator in equation Aw - f 
Cp pressure coefficient 
c coefficient in Trlcoml equation 
f right hand side in equation Aw •• f 


* This work was supported under NASA Contract ^ 1-14246, 
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F function in equation (5) 

G matrix iu equations (6), (7) 

I functional 

L Laplace operator 

Mdo freestream Hach number 

P pressure 

Q matrix operator 

R residual operator 

s streaoid.se coordinate 

t tlme-llke variable 

u,v velocity compo ne nt in x,y directions 

w dummy dependent variable 

x,y Independent variables 

a,6 coefficients in equation (1) 

Y ratio of specific heats .. 

5 Incremental difference operator (6a ” a**^ - a°) 
A finite difference operator (Ax ” x, - x^) 

e coefficient in equation (1) 

X parameter in equation (4) 

y switching factor, equation (34) 

V coe' 'Iclent in equation (1) 

p density 

p artificial density 

4> potential function 

* transpose operator 


TIME-DEPENDENT METHODS 


The Euler equations as well as Che unsteady full potential equation are 
hyperbolic in time, whether or not Che flow is locally subsonic or supersonic 
in space. With this in mind, Wellford and Hafez (ref. 1) solved the following 
augmented system for the transonic small-disturbance equation: 


ctu ■ ((1 - M^)u - m5u^) + V + v.u “ e.u 

t « 2 » x y Ixx 1 




u - V + v«v - e»v 
y X 2 yy 2 


( 1 ) 


where regularization terms 3v^ , v.u , v.v , e.u , e.v are added 

t’lxx2yy*12 

explicitly. The coefficients e, V, 3 are assumed to be small. An implicit 
Crank-Nicholson finite difference scheme was used for the time derivatives and 
a standard Galerkin finite element approach was used for Che space derivatives. 
Stability and convergence were rigorously analyzed. It was shown that a mini- 
mum amount of viscosity V^ , related to the magnitude of the drag, is needed. 

Results shown in fig. 1 are generally in agreement with finite difference 
results, but the accuracy is not acceptable for the element density used. The 
method converged slowly, but it may be useful for unsteady problems. 
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Hafez, Wellford and Murman (rof o\ jj 

to the full potential equation in thi cLe“«"e 


Pt “ -(P“)X - (Pv)y 


■ w 


where 


w - pY 1 _ (1 _ ^2 ^^2 ^^2 ^ 


Again, artificial vlscositv terms sc uah j 

added to the system. If we are damping terms (otw) are 

the transient behavior Is not Important and°a^f solution, 

possible* The work Is stlU InTSrew?”*^ ® Iterative procedure may 'be 

unsteady Euler**eJSatloni’^^ing^*‘Jlme^^^ element solutions of the 

non-conservative form with th! dependent variablH^bel' e^l^atlons are In 
elements are used along with a Galerkln Isoparametrlc 

result from their stud? Is sh^ ir^wf? S" derivatives. A 

but the rate of convergence Is results look quite good, 

~luiio„. Of the f”entill eq^tioS”'''”'’ 

VARIATIONAL METHODS 

many Impressive f Initrelem^t^aJlStJons^ elliptic and there are 

Bateman variational principle For fra«I ? 5,°^ *»ased on the 

variation ceases to L posftl^e 

Introduced mixed variational formulations Wellford and Murman (ref. 2) 

first given in terms of X u a^^i^is:"^ different functionals. The 


I(4>.u,v) 


where 


results shown in figure 3 for the cas^"*f**^ P ■ P(u,v) . Numerical 
Indicate approximately the right solution approximation 
second functional given In teras of^ J and rir”^ apparent inaccuracies. 


I$\ 


- u)2 + . ^)2i 


2 * 

(u - U(|)^) + (v^ - V(|)^) 




where 


-c. + 

1 Y-l 


' + r n a. _ _ . . _ 


The two associated Euler equations are the 


rrmh-f 4 *. j e^uacions are the 

hou„d:^'L“mSs ~:r:jreuLfL“i5? -a 

ditions across the shock (1 e mass^ls ^ admits the right Jump con- 

the definition of weak solu^l^i^s ?he tw^r^ ^ consequence of 

simultaneously. The method appears attwcti^e^^^ur^ "? 

able. As an extension to the above ideno ’ are yet avail- 

construct a general gradient mShoS ’ "" 
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( 6 ) 


0-0 

where G is 2 x 2 matrix. Notice if we choose G 
At 

*1 o and , 

i *2 ■ At , the system (equations 


(-■) 


(6) and (7)) 


(7) 

where 

is hyperbolic 


in time for subsonic as well as supersonic regions. Different choices of G may 
lead to faster convergence. Dissipation terms will be needed for stability and 
to handle flows with shocks. 


Eberle (ref. 4) has recently Introduced a variational finite element 
method for the full potential equation. Artificial viscosity terms are added 
which may be interpreted as an artificial compressibility. The density Is 
retarded upstream by a small (Ax) amount. The equations are solved in a 
transformed plane using a relaxation method. A calculated result is shown in 
figure 4. Other results in Eberle 's report indicate a noticeable dependence of 
the solution on element density, element shape and the magnitude of the 
artificial viscosity. We note that the work of Eberle motivated the artificial 
compressibility work discussed later. 


LEAST SQUARES AND OPTIMAL CONTROL METHODS 

Given the differential equation Aw * f, an associated functional whose 
second variation is positive definite can always be formulated using least 
squares, namely minimizing the residual 

||Aw - f||^ = (Aw,Aw) - 2(f,Aw) + (f,f) 

■ (A*Aw,w) - 2(A*f,w) + (f,f) (8) 

where the usual notation of inner products, norms, and adjoints are used. 

Notice that I - (A*Aw,w) - 2(A*f,w) is the Ritz variational functional for 
the problem A*Au * A*f which is automatically self-adjoint. However, the 
order of the equation has doubled, and therefore, a gradient method such as 

6w ■ - ■ - A*(Aw - f) (9) 

will converge very slowly. In order to accelerate the convergence of the 
iterative process modified gradient methods can be tised. We will discuss here 
three examples. 

In the first gradient method, a modification is made such that a Poisson's 
equation is solved at each iteration; namely 

L6w - I^ or 6w - l"^A*(Aw - f) . (10) 

Notice that higher order inter-element continuity of the shape functions is 
required for the least squares method compared to a Galerkin method. This 
problem may be avoided by writing the differential equation as a system of 
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lower order equations (see Lyim and Arya(refs. 5 and 6)). As an exasple, 
consider the Tricoml equation: 


Minimizing 
leads to: 


A(fr - c(x,y)<))^ + (j,^ - 0 


A.A* - 

2 

(c ^ ) + d) + (c(b ) + (c(b ) - 0 

^xx'xx ^yyyy ^ ^yy'xx '■ ^xx'yy '' 


Alternatively, 


and minimize 


let u ■ (|) , V ■ (j) ; hence cu + v ■■ 0 

X y X y 


+ (v - + Vy)^ I dA 


(11) 


( 12 ) 


(13) 


(lA) 


For simplicity, X is taken to be unity. In general, the choice of X may 
Improve the rate of convergence with respect to the mesh size as well as with 
respect to iteration. The Euler equations then may be written as 


(b + (b * u + V 
^xx ^yy X y 

2 

(c u ) + (cv ) - u ■ - d) 

XX ' y'x 


(cu ) +(v) - v"-d) 

X y y y ^y 


(15) 

(16) 
(17) 


A straightforward Iterative procedure to solve equations (15) to (17) is 
as follows: Given (j> , find u and v from equations (16) and (17), then 

using the most recent values of u and v in equation (15), a new value of (f) 
is obtained. This procedure can be described by 


(fiu) + (6v) 

A y 


(c u ) + (cv ) + (v ) + (cu ) 

x^xx y'xx ^ y^yy ^ x^xx 


(18) 


In terms of <|> , equation (18) is essentially (except for compatibility terms) 


L6(j) ■ A*A(|) + ••• (19) 

It should be mentioned that the advantage of using equation (14) Instead of 
equation (11) is that the elements used in the approximation of the latter must 
have continuous first derivatives, a fact which eliminates virtually all of the 
Important practical elements. It is also worth mentioning here that a linear 
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element approximation leads to undesirable convergence characteristics with 
respect to mesh size as discussed by Lynn and Arya (refs* 5 and 6)* It is 
obvious that if linear element trial functions are used for u and v , a 
bilinear approximation is needed for • 

Different least square formulations of the full potential equation are 
given in Table I. Chan and Brashears (ref. 7) solved the transonic small- 
disturbance equation by least squares similar to equation (11). Fix and 
Gunzburger (ref. 8) as well as Glowinskl et al. (refs. 9 and lu) used a formu- 
lation similar to equation (14). 

For the second method • we note that A is a second order operator. The 
"closest" positive operator to A*A is the Blharmonlc L*L » hence an iter- 
ative procedure based on 

L*L6w - - I - - A*(Aw - f) (20) 

w 

is expected to be fast provided the Inverse of (L*L) is easily obtained. 
Equation (20) is also equivalent to 

R - Aw - f , L*Z " A*R , L6w » Z . (21) 

The third method results from a more elegant decomposition of equation 
(20) obtained by modifying the inner product, namely 

I - ((Aw - f) , Q(Aw - f)) 

- (A*QAw,w) - 2(A*Qf,w) + (f,Qf) (22) 


hence g 

- A*QAw - A*Qf (23) 

where Q* ■ Q the choice of Q ■ makes A*L~^A effectively second order. 

The gradient method gives 

6w ■ A*L ^(Aw - f) (24) 

and the modified gradient method gives 

L6w - - A*L"^(Aw - f) (25) 

which may be decomposed into 

LZ - Aw - f , L6w “ -A*Z (26) 

Effectively, each step has an operator of zero order. Notice that 

I - - ((Lz) , L’^(LZ)) , - - (LZ,Z) . (27) 
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Upon Integration by parts, I becomes 

1 - (VZ,VZ) - ||vz||^ “ 

Recently Glowinski et al. (refs. 9 and 10) solved the full pov ntlal transonic 
flow equation using an equivalent procedure, namely 

minimize ^ w.r.t. w (29) 

under the constraint 


LZ “ Aw - f where Z ■ (j) - w . (30) 

The discussion with Professor Antony Jameson of Courant Institute helped the 
author in the preparation of this section of this pappr. 


TREATMENT OF SHOCKS 

In the above formulations, expansion as well as compression shocks are 
admitted since the potential flow Is reversible. In order to exclude expansion 
shocks, Glowinski et al. (refs. 9 and 10) introduced an entropy constraint* to 
the least squares formulation, namely 

V^(j) = M^(|) < + «• (31) 

s s 

A result shown in figure 5 indicates quite good agreement between the results 
and an exact solution for a shock~free airfoil « An alternative procedure has 
been Introduced by Bristeau (ref. 12) in which artificial viscosity terms are 
explicitly added to the continuity equation in addition to the above con- 
straint. A result is shown in figure 6. We notice here that the special 
assembly procedure used by Chan and Brashears (ref. 7) may effectively produce 
some sort of artificial viscosity in an obscure way without which the solution 
cannot be obtained. 


FINITE VOLUMES 

Finite volume methods use general nonorthogonal coordinates and consider 
the governing Integral equations as balances of mass , momentum, and energy 
fluxes for each finite volume defined by the intersection of the coordinate 
surfaces. Rizzl (ref. 13) applied the finite volume method to the Euler 
equation for transonic flows and calculated the time-accurate solution until it 
converged to a steady state. Factored explicit and implicit difference oper- 
ators were used to accelerate the calculations. Several computed examples are 
given in the reference. 


*A different procedure was used by Chattot (ref. 11) in his least squares for 
mulation of Euler Equations to exclude the expansion shock by imposing 

+ VSy > 0. 
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Jameson and Caughey (ref. 14) have recently applied the finite volume 
method to the steady full-potential equation in conservation form by using 
mixed-type flux operators. Centered difference operators are used with 

in supewonlc cells. Isoparametric mesh system is 
u.,ed and the equations are solved by relaxation. The method combines the 
advantage of finite elements for handling complicated geometries and the advan- 
tage of using simple finite difference schemes in the transformed coordinates. 

A typical calculation is shown in figure 7. 


The same concept has been recently used by Lucchi (ref. 
highi r order elements and with the velocities and the density 
used direct inversion to solve the matrix equation. 


15) with different 
as variables . He 


ARTIFICIAL COMPRESSIBILITY METHODS 

work of Eberle (ref. 4), we have explored the application of 
n artificial compressibility approach for both finite element and finite 
difference methods for mixed equations. The idea behind these new methods is 
o modify the density (and/or the speed of sound) in the supersonic region 
slightly (within the same order of the truncation error) and solve the 
resulting problem iteratively with standard methods used for the solution of 
elliptic problems . The density modification can be interpreted as an arti- 
ficial viscosity effect. The modified equation reads 


(P^^>x = (P4>y)y = 0 

For example, Jameson's fully conservative schemes 
form where 


(32) 

can be approximated by this 


P = P + y (uu Ax + w Ay) 

U = max (0, 1 - 1/M^) . 

We have tested an alternative form 


(33) 

(34) 


where 

Other possible, 
sound . 


P - yp As 
s 


Asp = ^ p Ax + - p Ay 
s q X q y ^ 


(35) 

(36) 


forms are under study including modification of the speed of 


Standard central difference formulas are used for equation (32) and S 
evaluated at the previous iteration. In the formula for p , upwind differ- 
encing is used in the supersonic region. Various Interatlve methods have been 
tested successfully including SOR (vertlcal-horizontal-SSOR) , /jl, a second 
order explicit method, fast solver and multigrid method. * ..amerlcal remit is 

I®'"®! explicit scheme. The rate of convergence 
<•■111 calculation is comparable to the relaxation methei. The method is 
still under development and details will be reported in a separate paper. 
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CONCLUDING REMARKS 

A considerable amount of work has been reported during the last few veers 
e]q>lorlng the application of finite element methods to transon'i<« 

In general, it is found that classical finitrli^L; problems, 

•ppllcbu. Modlflcelon. nst be JS. 

ge^Ji^ •ttrectlye <uid ere being eppllej to reelletlc elrfoil 
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TABLE I.- LEAST-SQUARE FORMULATIONS FOR FULL POTENTIAL EQUATION 

S - (pu)„ + (pv)„ 


n - Uy - 


>• ■ (p"-‘ - 


Functional 


I(u 


,v) - ffi 


where p - 1 “ ^ - 1)^”^ 


I(u 


,V,p) 


2 ^ o2 j. 1,2 

s + w + h 


where 


I(<t> 


,u,v) ”J‘J ‘ + (u - + (v - ((> 


y 

where 


,2 1 . ^ ^2 _ 


p - 1 - ^ M^(u^ + - 1)^"^ 


I((|>,u,v,p) - JJb^ + h^ + (u - <|)^)^ + (v - 4>y) 


(V 
where 


^ (u^ + - 1 ) 
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Figure 3.- Solution of transonic small-disturbance equation for flow past 
parabolic arc using mixed variational finite-element formulation of 
Hafez, Wellford, and Murman (fig. from ref. 2). 




t 



Figure 4.- Finite-element solution by Eberle for flow past NACA 0012 airfoil. 
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Figure 7.- Finite-volume solution of full-potential equation by Jameson and 

Caughey. 



Figure 8.- Solution of full-potential equation using artificial compressi- 
bility method. 



